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In this Supporting Information, we provide background information on our 
model of the in vitro Kai system and the calculations that we have performed. 
We will closely follow the outline of the main text. 

1 Allosteric Model 

In this section, we discuss in more detail the allosteric model of section I of the 
main text. We first present a statistical-mechanical description of the allosteric 
model. This allows us to describe the thermodynamics of the phosphorylation 
cycle. We then present a model based on concepts of transition state theory that 
allows us to describe the dynamics of the phosphorylation cycle, in particular 
the dynamics of the conformational transtions. Lastly, we briefly discuss how 
we have performed the simulations on the model in Fig. IB of the main text, the 
results of which are shown in Fig.lC of the main text. 
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1.1 Thermodynamics: A statistical- mechanical model 

The allosteric model relies on the following key assumptions: 

1. Each of the N = 6 monomers of a hexamer is either in an active or in an 

inactive conformational state. 

2. The conformations of the monomers are strongly coupled, such that all 
monomers of a hexamer are in the same conformational state at all times. 

3. Both in the active and inactive state, monomers can be (de)phosphorylated 

and (un)bind nucleotides. Wc assume that nucleotides have a higher affin- 
ity for monomers in the active state than for those in the inactive state. 
Consequently, nucleotide binding enhances the stability of the active state 
with respect to the inactive one. In contrast, we assume that phosphory- 
lation favors the inactive state. 

4. Nucleotide exchange is faster than phosphorylation and therefore in ther- 
modynamic equilibrium on the time scale of phosphorylation. 

The model makes the following further assumptions that are of secondary im- 



1. Each monomer has two phosphorylation states, phosphorylated and un- 
phosphorylated. Each hexamer thus has = 6 phosphorylation sites. 

2. Phosphorylation of the different monomers of a hexamer occurs sequen- 
tially around the hexamer. 

3. The unphosphorylated monomers can bind ATP, while the phosphorylated 
monomers can bind ADP. 

This leads to the following partition function for a hexamer that it is in 
a conformational state a, at phosphorylation level p, with q ATP and r ADP 
molecules bound: 



Here, /3 is the inverse temperature, is the energy of an unphosphorylated 

monomer in state a (with no nucleotide bound), and K^''^ is the dissociation 
constant for the binding of species s to the hexamer in state a; T denotes ATP, 
D ADP, and P a phosphate group Pi. Since nucleotide exchange is assumed to 
be fast, it is meaningful to integrate over the number of nucleotides. This yields 
the following partition function for a hexamer in state a with phosphorylation 
level p: 



portance: 




(1) 
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Figure 1: The chemical potential of a KaiC hexamer as a function of the phos- 
phorylation level p, for both the active and inactive conformational state, and 
for the symmetric model of Fig.lA in the main text. The chemical potential 
fj," (p) is given by Eq. [3l 



The (excess) chemical potential of a hexamer in conformational state a at 
phosphorylation level p is given by 

fi^ip) = -kBT\n[Z^{p)]. (3) 

We stress that the chemical potential of a KaiC hexamer depends upon the 
chemical potentials (the concentrations) of the nucleotides: ^"(p) = fJ.°'{p; fJ.^ , fjP,^^). 

We will now consider the symmetric model of Fig.lA of the main text. In 
this model, the energy levels of the active and inactive conformational state are 
mirror images of each other. This is for reasons of clarity, and not because it 
is essential. In fact, the full model of section III is asymmetric, with the active 
state being more stable than the inactive one. 

Fig. [1] shows the chemical potentials of a KaiC hexamer in the active and 
inactive state, respectively, as a function of the phosphorylation level, for the 
energy diagram shown in Fig.lA of the main text. Here, fJ."{p = 6) — /i"(p = 0) 
corresponds to the free-energy change upon fully phosphorylating a KaiC hex- 
amer in conformational state a at constant chemical potentials of the ATP, ADP 
and Pi molecules, but without an ATP hydrolysis reaction (we thus consider 
the reaction KaiC + Pi ^ KaiCPi) . It is seen that the free energy increases 
markedly for both conformational states, meaning that the probability that a 
hexamer would fully phosphorylate spontaneously, is essentially zero. Indeed, 
the essence of our allosteric model is that in the active conformational state, the 
energy from ATP hydrolysis is used to phosphorylate the KaiC hexamer, while 
in the inactive state dephosphorylation occurs spontaneously. 
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When ATP is hydrolyzed to p-fold phosphorylate a single KaiC hexamer in 
the active state, the total change in free energy of the system is 

AG^ip) = fi^ip) - M^(0) +p(a*° + - /i^), (4) 

where H^ip) is given by Eq. [31 If in the active state binding of both ADP 
and ATP is strong (i. 
expression reduces to 



and ATP is strong (i.e. [ADP]/Kd'^°^, [ATP]/Kd''^'^^ > 1), then the above 



AG^{p) = p(-AG^'T + AG^'P + AG^'D + AGhydro) 

— -P^G'^;AATP^ApADP- (5) 

Here, AG"''' — +kBT\nK'^''' is the binding free energy of species s and AGhydro 
is the standard reaction free energy of an ATP hydrolysis reaction. The overall 
free energy change AG^(p) corresponds to that of p phosphotransfer reactions 
AATP ApADP on the active KaiC hexamer (A denotes a subunit in the 
active state). The free-energy change be understood by noting that in the limit 
of strong nucleotide binding considered here, the unphosphorylated monomers 
are essentially always occupied by ATP, while the phosphorylated monomers 
are essentially always occupied by ADP; see also Fig.lA of the main text. 

Fig. [5] shows the free energy of the system in the presence of ATP hydroly- 
sis. In the active conformational state, KaiC binds ATP. ATP hydrolysis then 
drives the phosphorylation of the hexamer, and the reduction in free energy 
of the whole system is given by Eq. [SJ When the hexamer is (nearly) fully 
phosphorylated, it flips to the inactive conformational state. In the inactive 
state, nucleotide binding is weak and, as a result, ADP is released. The hex- 
amer now dephosphorylates spontaneously; since nucleotide binding is weak, 

the free-energy change is given by B/cbT In [PiJ/ifp^ (see Eqs. [5] and [3]). At 

low phosphorylation levels, the inactive hexamer flips back towards the active 
conformational state. KaiC rebinds ATP and the phosphorylation cycle starts 
over again. After one full phosphorylation cycle, the free energy of the system 
has been reduced by the free-energy change corresponding to 6 ATP hydrolysis 
reactions: AG = 6 (AGhydro + fceTln ([ADP][Pi]/[ATP])). 

1.2 Dynamics: a transition-state theory of the conforma- 
tional transitions 

So far, we have discussed the thermodynamics of the phopshorylation cycle. We 
will now discuss the dynamics of the cycle, in particular the dynamics of the 
conformational transitions. This is important, because while the large amplitude 
oscillations as observed experimentally require that the hexamers should not 
flip at intermediate phopshorylation levels, the stability of one conformational 
state with respect to that of the other, does change in sign at intermediate 
phosphorylation levels: in the symmetric model considered here and in Fig.lA 
of the main text, the active state is more stable for p < 3, while the inactive state 
is more stable for p > 3. How can we explain that the conformational transitions 
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Figure 2: The free energy of the system as a function of the phosphorylation 
level in the presence of ATP hydrolysis for the symmetric model of Fig.lA in 
the main text. The solid lines denote the path of the system. Driven by ATP 
hydrolysis, a KaiC hexamer is phosphorylated in the active state. When the 
KaiC hexamer is (nearly) fully phosphorylated, it flips from the active to the 
inactive conformational state. In the inactive state, ADP is released and the 
hexamer dephosphorylates spontaneously. At low phosphorylation levels, the 
hexamer flips back to the active state. The active hexamer rebinds ATP and 
the phosphorylation cycle starts over again. The dotted lines correspond to 
driven phosphorylation of inactive KaiC and spontaneous dephosphorylation of 
active KaiC. 

predominantly occur when the hexamers are either nearly fully phosphorylated 
or fully unphosphorylated? 

This is a difficult question to answer, because it requires knowledge of the 
microscopic dynamics of the transition paths between the conformational states. 
However, if we assume that nucleotide binding is an important component of 
the reaction coordinate that describes the conformational transitions, then we 
can make an estimate of the flipping rates using a mesoscopic model based on 
concepts from transition-state theory [1]. 

If nucleotide binding contributes to the reaction coordinate of the confor- 
mational transitions, then we cannot integrate it out as we have done so far. 
To derive the flipping rates, we start by considering the free-energy difference 
between two conformational states with the same number of nucleotides bound: 



Here, A£^,„ = E^ — E]^ and lS.Es — k^T In K^"'' / . If we assume, for simplic- 
ity, that the model is symmetric, Aii^„i = 0, and that the difference in binding 
energy between the active and inactive state is the same for ATP and ADP, 



AG(p, r) = N/SE^ + pAEp + qAEr + rAEu ■ 



(6) 
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Ai?T = AE'd = Ai?T/Di then the above expression reduces to: 

AG{p, n) = pAEp + hAEt/d, (7) 

where n is the number of nucleotides that are bound. We iterate that phospho- 
rylation favors the inactive state, and hence AEp > 0, while nucleotide binding 
favors the active state, Ai5x/D < 0. We can use the above expression to es- 
timate the flipping rate if we assume that nucleotide binding is the dominant 
reaction coordinate for the flipping process. 

This is illustrated in Figs.[n]and[31 for the case \AEp\ = -|Ai?„|. Fig.[3]shows 
a sketch of the free-energy surface AG{p = 3, n, c) of a three-fold phosphorylated 
hexamer as a function of the number of bound nucleotides, n, and as a function 
of an order parameter that describes the conformational state of the hexamer, 
c; the parameter c is zero if the hexamer is in the inactive state and one if it 
is in the active state. Clearly, we do not know what would be the best order 
parameter to describe the conformational transition, let alone what the free 
energy would be as a function of this order parameter for different values of 
n. Nevertheless, we do have some knowledge of the free-energy surface: we 
know how the free energy AG{p, n, c) changes as a function of n for c = 
and for c = 1 - this is given by the free energy of nucleotide binding to the 
inactive and active state, respectively; this free energy is related to the log of 
the partition function in Eq. [TJ We therefore make the minimal assumption that 
the free-energy surface AG{p, n, c) is given by a linear interpolation between the 
two functions AG{p,n,0) and AG{p,n, 1). This leads to the surface shown in 
Fig.l 

In the active state, nucleotides bind the hexamer very strongly and, conse- 
quently, n « 6, while in the inactive state nucleotides bind the hexamer rather 
weakly and n « 0. The two (meta)stable states of the hexamer are thus the ac- 
tive state with six nucleotides bound and the inactive state with no nucleotides 
bound. These two states are separated by a "transition-state" surface: in or- 
der to go from one (meta) stable state to the other, the system has to cross a 
free-energy barrier. We assume that the transition state is given by the saddle- 
point in the free-energy surface AG{p,n,c), as shown in Fig. [31 This means 
that both the location and the height of the free-energy barrier for flipping are 
determined by that number n* for which the two states become equally stable, 
AG{p,n*) — (see Fig. |3] and Eq. [7]). Clearly, the location of the transition 
state depends upon the phosphorylation level p of the hexamer: in the symmet- 
ric model considered here, the two conformational states are equally stable if 
the number of bound nucleotides is n* = p (see Fig. |4]and Eq. [7]). The height 
of flipping from the active to inactive state is thus given by 



/3AGl^i(p) = -ln 



Z^iP. r)Siq + r - p)/Z^{p, 0, 6) 



(8) 
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Figure 3: The free energy of a KaiC hexamer for a phosphorylation level of 
p = 3, as a function of the number of bound nucleotides n and as a function 
of an order parameter c that denotes the conformational state of the KaiC 
hexamer: it is zero if the hexamer is in the inactive state and one if it is in 
the active state. In the inactive state, essentially no nucleotides are bound, and 
n « 0, while in the actives state, because of strong nucleotide binding, n « 6. It 
is seen that in order to flip from the inactive to active state, the system has to 
cross a free-energy barrier; the transition state is denoted by the red cross. We 
imagine that the height of the barrier to go from the inactive to active state is 
given by the free-energy to add three nucleotides, while the barrier to flip from 
active to inactive is given by the free energy to remove 3 nucleotides. 

while the barrier height for the reverse transition is given by 



Here, Z°'{p, q, r) is given by Eq.[TJ In words, if a p-fold phosphorylated hexamer 
is in the active state, with 6 nucleotides bound, then in order to flip to the 
inactive state with no nucleotides bound, it has to cross a barrier with a height 
that corresponds to the energetic cost of removing 6 — p nucleotides. Conversely, 
the height of the barrier for an inactive hexamer, with no nucleotides bound, 
to flip to the active state, is given by the free energy to add p nucleotides. 
Neglecting entropic factors, the height of the free-energy barrier thus scales 
linearly with the phosphorylation level, leading to the exponential flipping rates 
of Eqs.l and 2 of the main text. 

1.3 Numerical calculations on the allosteric model 

The chemical reactions of the model in Fig. IB of the main text are: 



/3AGi*_^a(p) = - In ^ Z\p, q, r)Siq + r - p)/Z\p, 0, 0) 



(9) 
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Here, Ct corresponds to an active KaiC hexamer with phosphorylation level 
i, while d corresponds to an inactive KaiC hexamer with phosphorylation level 
i. The first, reversible reaction corresponds to the conformational transitions of 
the KaiC hexamers with forward and backward flipping rates fi and bi, respec- 
tively, the second corresponds to phosphorylation of active KaiC at rate kps, 
while the third reaction corresponds to dephosphorylation of inactive KaiC at 
rate /cdpg. 

To study the phosphorylation behavior of a single KaiC hexamer, we cannot 
use macroscopic rate equations based on the law of mass action: these equations 
would correspond to the average of a population of KaiC hexamers. To simulate 
the behavior of a hexamer, we have performed kinetic Monte Carlo simulations of 
the zero-dimensional chemical master equation corresponding to the reactions in 
Eq.[T0]|2|. The solid line in Fig.lC corresponds to the results of those stochastic 
simulations. 

To study the time evolution of the average phosphorylation level of an en- 
semble of KaiC hexamers, we have used macroscopic rate equations based on 
the law of mass action. The chemical rate equations that correspond to Eq. [TOl 
are: 

kp,[a_l] - fcp,[a] + 6,[C,;] - MC] (11) 

hps [a-i] - hps [a] + n [a] - h [q] (12) 

The dashed line in Fig.lC of the main text corresponds to the numerical results 
of propagating these ordinary differential equations. 

The results in Fig.lC of the main text were obtained with the following 
values for the parameters: fcps = O.Olhr""'^, fcdps = 0.05hr~^, fi — 0.1^~* hr~^ 
and hi = O.V hr"^. We note here that the above chemical rate equations, 
corresponding to the model in Fig. IB of the main text, are identical to those 
of the simple model with differential affinity as described in section II with 
[A]t = (see Eqs. [TTIIT^ . Also the values for the parameters are indentical to 
those in the differential-affinity model of section II (see section below). 

2 Simple Models with Differential Affinity 

In this section, we first provide background information on the simplified model 
of the Kai system discussed in section II of the main text. We then briefiy 
discuss a more generic class of differential affinity models. 



4Q] 

dt 

d[C^] 
dt 



8 



Figure 4: The free energy of the active and inactive state as a function of 
the number of bound nucleotides, n, for three different phosphorylation levels, 
p = 0,3, 6. In the active state, n « 6, while in the inactive state, n 0. The free 
energy of p-fold phosphorylated hexamer in state a with n nucleotides bound is 
given by AG"(p, n) = — fceTln^^ ^ ^"(p, 9, r)6{q + r — n), where Z"{p, q, r) is 
given by Eq. [1] 

2.1 A minimal differential affinity model of the Kai system 

In the simple differential affinity model of section II of the main text, we assume 
that only a single KaiA dimer can bind to an active KaiC hexamer. The chemical 
reactions of this model are given in Eqs.3-5 of the main text. They correspond 
to the following mass-action kinetic equations for the concentrations [Cj] and 
[Ci] of KaiC in the active and inactive states, [AC;] of the KaiA-KaiC complex, 

fcpf [AC,.i] - fc^f [A][C,] + k^''[AC,] 
+6^Mbo[Go]-S,,ef6[C6] (13) 
-kpf[AC,] + k^'[A][C,]~k^''[AC,] (z^6) (14) 

fcdps([C,-i] - [C,]) - S^^oboiCo] + S^,ef6[Ce] (15) 

-[A] ^ kf[G,] + ^(fcf + A:pf )[AC.] , (16) 

4=0 i=0 

where 5ij is the Kronecker delta, fcdps is the spontaneous dephosphorylation 
rate, fcpf is the rate of phosphorylation catalyzed by KaiA, and fi and hi are 
the flipping rates as defined in Eqs.l and 2 of the main text. The rates of KaiA 
binding to and unbinding from active KaiC are respectively k^^ and kf^^, with 



and [A] of free KaiA: 

dt 

d[AC,] 
dt 

d[C,] 
dt 

4M 

dt 
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the latter dependent on the number of i of phosphorylated monomers in the 
KaiC hexamer. Because we generally choose parameters such that (un)binding 
of KaiA to KaiC is much faster than (de)phosphorylation, it is an excellent 
approximation to assume that these binding reactions are equilibrated. In this 
case, we do not have to keep track of [Ci] and [AC^] separately; instead, we 
obtain dynamical equations for the total concentration [C^Jt = [Ci] + [ACi] of 
KaiC in the active state and for [C^]: 



d[C,]T fcpf[A] fcp/[A] 

- [Ci_ljT ~ , r . 1 [C. 



JT 



dt K,^i + [M K, + [AY 

+S^,ob^[C,]~S,^ef^[C^]T (17) 



dt 



fcdps([C,-l] - [C,]) - SMC,] + (5,,6/,;[C,]t (18) 



along with a constraint equation giving the free KaiA concentration [A] implic- 
itly: 



Here [A]t is the total concentration of KaiA and Ki — kf^ /k^^ is the disso- 
ciation constant for KaiA binding to i-fold phosphorylated KaiC in the active 
state. We implement differential affinity by setting Ki = /3a' with a > 1. Here, 
and in the next section, the differential equations were solved using Matlab. 
We note that the assumption that KaiA binding is fast is convenient for some 
purposes but not essential; we have verified that the model's predictions are the 
same whether we use Eqs. [T51-[TBlor the reduced set Eqs. [171 [T^ 

PARAMETER VALUES GO HERE (or probably, better yet, in a table). 



2.2 Parameter dependence and bifurcation behavior 

It is natural to ask how the model's behavior changes as the parameters are 
varied from the specific values listed in the previous subsection**OR: IN THE 
TABLE**. This task is facilitated by its simple cyclic structure, which allows 
one to prove that the system of equations 1 1 3I[TB1 alwavs admits exactly one fixed 
point. As the total concentration [A]t of KaiA is increased from zero, this fixed 
point becomes unstable through a supercritical Hopf bifurcation. The resulting 
stable limit cycle persists through a fairly broad range of [A]t values before 
finally disappearing at a second supercritical Hopf bifurcation where the unique 
fixed point regains its stability. This behavior can be understood based on argu- 
ments similar to those put forth in Section IIIC of the main text: If [A]x is too 
large, then the concentration of KaiA is no longer limiting, and differential affin- 
ity cannot act to synchronize oscillations. On the other hand, if [A]t beocmes 
too small (while the other parameters are held fixed), then the phosphorylation 
reactions will proceed too slowly compared to the dephosphorylation reactions. 
The first phosporylated hexamers will then be dephosphorylated and return to 
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state Co before the remaining KaiC complexes reach state Cg; the hexamers 
in state Co will win the competition for the limited KaiA, preventing the syn- 
chronized release of the remainder of the KaiC. Varying other parameters has 
similar effects. For example, as the parameter in Ki = f3a^ is increased, the 
limit cycle will eventually collapse into the fixed point in a Hopf bifurcation. 
Indeed, if the Ki become too large, then KaiC hexamers cannot efficiently bind 
all of the available KaiA, and differential affinity is no longer possible. Figure ?? 
shows a 2-parameter bifurcation diagram as a function of [A]t and f3. We will 
see that similar bifurcation behavior reappears in more realistic models of the 
Kai system. 

2.3 Generic differential ctfRnity model 

The simple model discussed above and in Section II of the main text can be seen 
as one example of a more generic class of models that use the same mechanism 
to synchronize oscillations. Here, we briefly discuss this broader perspective on 
differential affinity. A fuller mathematical analysis will appear in a forthcoming 
publication. 

We begin by considering the following cycle: 

C'^'O ki ^ ^ kN ^ /'on\ 

— > L^l — > . . . — > <^JV-1 — * ^JV — * ^0, (,^Uj 

where C; is a protein that has been i-fold covalcntly modified. At least two of 
the reactions require a catalyst A acting with Michaelis-Menten kinetics, and 
those steps that are not catalyzed by A are simple first-order reactions. 

Suppose that the reactions Co ^ Ci ^ . . . ^ Cj require the catalyst A 
(with j < N). Then, one might imagine that this system will oscillate if a) 
the concentration of A is sufficiently low and b) if the dissociation constants Ki 
for the binding of A to Ci satisfy Kq < Ki < . . . < Kj. These two conditions 
together ensure that A first binds to Co and catalyzes the reaction Co Ci; 
only when the concentration of Cq has dropped almost to zero does A begin to 
bind to Ci and catalyze the next reaction, and so on until state Cj is reached. It 
turns out, however, that these requirements alone are not sufficient. In addition, 
we must demand that c) the distribution of arrival times of different C molecules 
at Cj is not too wide compared to the average time to travel from Cj back to 
Co and d) that the distribution of arrival times back at Co is not too broad. If 
the former condition, c), does not hold, then the fastest C molecules will reach 
Co while some of the slower C molecules have still not reached C j . As we noted 
in the previous subsection, oscillations cannot survive such a situation: Because 
A binds most strongly to Co, the fastest C molecules will siphon away A from 
the slowest that still need A to progress around the cycle; these will then slow 
down further until all sychronizing effect of differential affinity has been lost. 
Similarly, differential affinity fails when the arrival of C molecules at state Co 
is too spread out. Then, there are relatively few C's competing for A molecules 
at any given instant, and even those with the highest number i of modifications 
can continue upwards towards Cj. Thus, for oscillations to occur, N — j must 
be neither too small nor too large. 
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We have succeeded in finding limit cycles in this model with j as small as 1 
(corresponding to A binding to the two states Co and Ci) and N as small as 
4. Fig ?? shows oscillations in such a model with **PARAMETER VALUES 
HERE**. Although the oscillations persist only up to TV = 7 for these same 
parameter values, there appears to be no maximum allowed value of N if the 
rate fcps of the transitions Cj Cj+i ^ . . . ^ C^v ^ Co is allowed to decrease 
as l/{N—j). Indeed, in this case, the mean time to travel from Gj to Cq remains 
constant, while the distribution of travel times becomes narrower, which only 
enhances oscillations. 

One could imagine many variations on the model just described. For exam- 
ple, we anticipate that for appropriate parameter values, oscillations will also 
occur when A binds to and is sequestered by Co but is not required for the 
transition Co — > Ci, while still catalyzing the reaction Ci — *■ C2. One might 
also consider cases in which differential affinity acts on two separated blocks of 
reactions Co ^ Ci ^ . . . ^ Cj and Cn — + C„+i — ^ . . . ~> Cn+k- In each case, 
the same basic principles of differential affinity should be at work. 

The models of the Kai system in the main text are related to this generic 
class of models. The phosphorylation cycles of both the minimal differential 
affinity model of section II and of the full model of the Kai system in section III 
have an active branch where KaiC is phosphorylated, and an inactive branch 
where KaiC is dephosphorylated. Also in both models, KaiA catalyzes the 
phosphorylation reactions on the active branch. Both ingredients are inspired 
by experimental observations; together they give a concrete example of how 
the abstract cycle of Eq. [20] might be implemented. However, the discussion of 
the generic model above shows that from the perspective of synchronising the 
oscillations, there is no need to make a distrinction between an active and an 
inactive branch. Indeed, the same formalism could be applied to cycles made of 
more than two allosteric conformations, or of multiple different sorts of covalent 
modifications, or created in any number of other ways. The differential affinity 
mechanism thus has the potential to be generalized far beyond the Kai system. 

3 Full Model of the Kai System 
3.1 Model: the chemical rate equations 

The chemical reactions that describe the full model of the Kai system are given 
by Eqs.7-11 of the main text. Using the law of mass action, this leads to the 
following set of macroscopic chemical rate equations for the concentration of 
[C], [AC,], [C,], [B™C,], and [A„B„,C,]: 
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fcps[C.,-i] + /,:dps[Ci+i]-(fcps + fcdps)[Ci] + A:pf[AQ_i] 

-MC^]+h[C^] - kf[A][C^] + fcf^[ACi] (21) 

fcf [A] [Ci] - kt\ACi] - fcpf [ACi_i] (22) 

fcps [Ci-i] + fcdps [Cn-i] - (fcps + fcdps) [Ci] + fi [Ci] - hi [Ci] 
-fcff[B]™[Q] + fcfb[B„Q] (23) 

fcps [BmC-i— i] + A^dps [BjTiCi-i-i] (A)ps + A^dps) [BrnCj] 
-^f [A]"[B„Q] + ~kt''[Am^mCi\ (24) 

A^ps [AmBj^Cj—i] +A;dps[AmBTOCj-|-i] — (fcps+fcdps)[ATOBj„Ci] 
[Al^p^Ci] - fef'^fA^B^Ci] (25) 
Here, the concentrations of free KaiA and KaiB, [A] and [B], are given by: 

N 

[A] = [A]T-^([AQ]+m[A„B„Ci]) (26) 

i=0 

AT 

[B] = [B]T-^(m[B„Q]+m[A„B„Q]) (27) 

i=0 

The phosphorylation and dephosphorylation rates on the active branch are fcps 
and fcdps, respectively, and the flipping rates are fi and hi. The active state 
can bind KaiA with forward and backward rates fc^' and kf'^ , and KaiA can 
catalyze phosphorylation with the rate A'pf. We assume that in the inactive 
state KaiC can bind m = 2 KaiB molecules with forward and backward rates 
kf^ and kf^, respectively. This KaiB-KaiC complex can then sequester m = 2 
KaiA molecules with forward and backward rates kf^^ and kf^. 

3.2 Model: The free-energy difference between the active 
and inactive state of KaiC 

In the presence of only KaiA, KaiC is phosphorylated to a very high level of 
90 — 95%. This requires that the active state of KaiC be more stable than the 
inactive one. However, in order for KaiC to sustain oscillations in the presence 
of both KaiA and KaiB, KaiC should be able to flip from the active state to 
the inactive one at higher phosphorylation levels. This means that the free- 
energy difference AG(p) between KaiC in its active and inactive conformational 



d[Ci 
dt 



d[ACi 
dt 

d[Ci 
dt 



d[BmCi 

dt 



d[AfyiBrnCi 

di 
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^pS; ^ps 


0.025 hr"^ 


^dps; ^dps 


0.4 hr"^ 


kpf 

fi 


1.0 hr"^ 

{10-^ 10-^ 10-^, 10-3, 10-2, 10-1, 10} hr-^ 


hi 


100 hr-^ 




1.72-10^ M-^hr-^ 




{10, 30, 90, 270, 810, 2430, 7290} hr"^ 


~kr 


2.97-1012 X {0.01, 1, 1, 1, 1, 1, 1} M-^hr"^ 




M02 X {10, l,l,l,l,l,l}hr-i 
2.97-101^ X {0, 1, 100, 100, 1,0,0} M-^hr-^ 




^Ab 


100 hr-^ 


[A]t 


0.58 ijM 


[B]t 


1.75 nM 


[C]t 


0.58 ijM 



Table 1: List of parameter values for the full model in Section III of the main 
text. 

states should be strongly negative, but also that it should rapidly decrease 
in magnitude as KaiC becomes fully phosphorylated. A simple model that 
captures this is one in which the addition of each phosphate group decreases 
the free-energy difference by an amount Ai?p, and moreover where the creation 
of an interface between a phosphorylated and an unphosphorylated unit also 
destabilizes the inactive state by an additional amount e: upon the addition of 
the last phosphate group, two such interfaces are removed, leading to an extra 
change in the free-energy difference of 2e in favor of the inactive state. This leads 
to the following expression for the free-energy difference between the active and 
inactive state: 

AG(p, q, r) = NAE^ + pAEp -e^ninj+ qAEt + rAEr, . (28) 

Here, rii denotes the phosphorylation state of unit i - rij = 1 if unit i con- 
tains a phosphate group and zero otherwise - and the sum includes all 
nearest neighbor units of KaiC. For the full model, AE„^, AE^^, AE^j < 0, and 
AEp, e > 0. For simplicity, we assume that AEp = 2e. This yields the follow- 
ing expressions for the transitions between the active and inactive states: for 
the forward rate we have fi = 5j^~'^ for 1 < i < TV and fi = 10 • 5j'^~'^ for 
i = 0,N (this is due to, respectively, the creation or removal of the interfaces 
between phosphorylated and unphosphorylated units). The backward rate bi is 
independent of i and much larger than the forward rate fi, so that the stability 
of the inactive state is only due to binding of KaiB. 
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3.3 Setting the parameters 

Using the expressions for fi and bi just derived, the full model contains 39 
parameters. However, their values are very much constrained by the large body 
of experimental data on this system. We now describe how we have determined 
the parameters, and how critical their precise values are for the behavior of the 
model. Unless indicated otherwise, the exact values of the parameters for the 
full model are summarized in Table [TJ 

Concentrations The concentrations of KaiA and KaiB dimers are [A]t = 0.58^M 
and [B]t — 1.75/j,M, respectively, and the concentration of KaiC hexamers is 
[C]t = 0.58/zM. This corresponds to a concentration ratio of (KaiA dimers) : (kaiB 
dimers):(KaiC hexamers) = 1:3:1. The corresponding monomer concentrations 
are 1.17/LtM KaiA, 3.5/LtM KaiB and 3.5/j,M KaiC. Oscillations in phosphoryla- 
tion have been observed for these concentrations in the in vitro experiments of 
Tomita et al. [3]. However, it should be noted that for the results in this article 
the ratios of KaiA and KaiB to KaiC, [A]t/[C]t and [B]t/[C]t respectively, are 
more important than absolute concentrations. Hence, we often express concen- 
trations in units relative to [C]t. 

Flipping rates Following the discussion in the last paragraph of the previous 
section on the model for the free-energy difference between the active and inac- 
tive state of the KaiC hexamers, we next discuss the values of the flipping rates. 
As discussed in that paragraph, it is important that: 1) in the absence of KaiB, 
the active state has a lower free energy than the inactive one; 2) the hexamers 
should not flip at intermediate phosphorylation levels. In the previous section, 
we presented a model for the free-energy difference between the active and in- 
active state. Here we give the values of the flipping rates that are consistent 
with this free-energy difference, which is requirement 1), and with requirement 
2). The forward rate is /, = 6-/^-' ior 1 < i < N and f^ = 10■ ^7^-* for 
i — 0, N, with 7 = 0.1 and S = 2hr^^. The backward rate bi is independent of 
i and given by 6 = lOOhr"^. As long as the two requirements of this paragraph 
are satisfied, however, the precise values of the flipping rates are not important 
for the behavior of the model. 

KaiC alone and KaiC + KaiB: Spontaneous (de)phosphorylation rates 

The spontaneous phosphorylation and dephoshorylation rates are chosen such 
that the phosphorylation behavior of KaiC in the absence of KaiA and KaiB 
and that in the presence of only KaiB agrees well with experiment (see Fig. 4 of 
the main text) . This yields the following rates for, respectively, the spontaneous 
phosphorylation and spontaneous dephosphorylation reactions, both for active 
and inactive KaiC: fcps — kps = 0.025hr~^ and fcdps = ^dps = 0.4hr~^. The 
identical rates for active and inactive KaiC ensure that the phosphorylation 
behavior of KaiC in the absence of KaiA and KaiB is the same as that of KaiC 
in the presence of KaiB (which stabilizes the inactive branch). The values of 
these rate constants are not free and have to be carefully chosen, because their 
sum determines the relaxation rate of the phosphorylation level of KaiC (in the 
absence of KaiA) , while their ratio determines the steady-state phosphorylation 
level of KaiC (in the absence of KaiA) . 
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KaiC + KaiA: rates of KaiA-catalyzed phosphorylation reactions The 

phosphorylation rate of KaiC in the presence of KaiA is determined by two 
factors: 1) the binding affinity of KaiA for KaiC; 2) the rate fcpf for the KaiA- 
catalyzed phosphorylation reaction. As discussed in the main text, both the 
mechanism of differential affinity and temperature compensation require that 
the binding affinities of KaiA for KaiC be high. Given these high binding affini- 
ties, the rate at which KaiC is phosphorylatcd in the presence of KaiA is mostly 
determined by kpf. This rate is thus chosen such that the phosphorylation rate 
of KaiC in the presence of KaiA agrees with experiment. This gives fcpf = lhr~^. 
This rate constant cannot be freely chosen, because it directly affects the phos- 
phorylation rate of KaiC in the presence of KaiA. 

The mechanism of differential affinity requires that the affinity of KaiA for 
KaiC be high, but decrease substantially as the phosphorylation level of KaiC 
increases. Provided that these constraints arc satisfied, the precise values are of 
less importance. We have chosen the following expressions for the forward rate 
fcf f and backward rate kf°: kf = 1.72 • lO'^M-ihr"^ and fcA^ = /3Q;*hr"\ with 
a = 3 and /? = 10. 

The binding of KaiA and KaiB to inactive KaiC The stoichiometrics of 

the complexes of KaiC bound to KaiB and KaiA are not known. We assume 
that inactive KaiC, Cj binds m = 2 KaiB dimers, while the complex B2Ci can 
sequester m = 2 KaiA dimers. The value of m = 2 is not very critical, as 
long as enough KaiA can be sequestered by inactive KaiC the amount of KaiC 
available to catalyze phosphorylation is limited. We find that for to > 2 the 
system also oscillates, although the phase boundaries, as shown in Fig. 6 of the 
text, do shift. Clearly, more experiments are needed to resolve the compositions 
of these complexes. 

Temperature compensation requires that the binding affinities of KaiB for 
KaiC be high. But as long as this requirement is fulfilled, the precise values 
are of less importance. We assume that KaiB binds inactive KaiC with forward 
rate = 2.97- lO^^ ^ {o.Ol, 1, 1, 1, 1, 1, IjM-^hr"^ and backward rate kf^ = 
1 • 10^ X {10, 1, 1, 1, 1, 1, l}hr-i. The low affinity of KaiB for Co is to ensure 
that unphosphorylated hexamers that are in the inactive state can rapidly flip 
towards the active state. 

As explained in the main text, it is important that the affinity of KaiA 
for KaiB is low for KaiB associated with the highly phosphorylatcd inactive 
KaiC hexamers, but increases strongly for KaiB that is bound to the Icss- 
phosphorylatcd inactive KaiC hexamers: the precise values of the dissocia- 
tion constants arc not critical. We have chosen the following forward rate kf^^ 
and backward rate fcf^: kf = 2.97 • 10^^ x {0, 1, 100, 100, 1, 0, OjM-^hr"^ and 
k^^ = lOOhr"^ The low affinity of KaiA for B2C0 is to ensure that unphos- 
phorylated hexamers that are in the inactive state can rapidly flip towards the 
active state. 

Varying the parameters with temperature In order to calculate how the 
oscillation period in our model varies with temperature, we would have to know 
how the rate constants would vary with temperature. This would require knowl- 
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edge of the activation barrier for each of the reactions, which we do not have. 
We can, however, test how sensitive the oscillation period is to changes in the 
rate constants. Here, one could choose various strategies: one could vary each 
of the rate constants individually, or one could vary all the rate constants simul- 
taneously, either in a correlated or an uncorrelated manner. We have performed 
many such tests, and they all reveal that the model is very robust to changes in 
temperature (given the experimentally observed insensitivity of the phosphory- 
lation rates to changes in temperature (see main text)). Since the dissociation 
constants and the flip rates enter the model in very different ways, we show in 
the main text how the oscillations change when we vary these two groups of 
parameters. 

When we change the dissociation constants of KaiA and KaiB binding to 
examine the effects of temperature changes, we simultaneously change the dis- 
sociation constants for KaiA binding, Ki and Ki, hy a. factor C, and the as- 
sociation and dissociation rates for KaiB binding, kf^ and kf^^ by l/^fC and 
VC, respectively. Because KaiB (un) binding is fast, only the ratios of their 
rates — the dissociation constants — matter, and this is a choice that is consistent 
with changing the dissociation constants by a factor C . 

3.4 Reduced model 

When KaiA binding and unbinding is sufficiently rapid, it is possible to further 
simplify Eas. [2Tll25l In this case, we only explicitly take into account the binding 
of KaiB to the inactive KaiC hexamers and we assume that both KaiA associa- 
tion to the active KaiC hexamers and KaiA sequestration by KaiB bound to the 
inactive hexamers is very rapid and can be treated as in chemical equilibrium. 
This leads to the following reduced set of macroscopic chemical rate equations 
for the total concentration [Cijx of KaiC in the active state, both with and 
without bound KaiA, the concnetration of KaiC in the inactive state, [C J , and 
the total concentration [BmC^Jx of KaiC in the inactive state with KaiB bound, 
both with and without bound KaiA: 



d[C,]T 



arjC._i]T+4^1C^i]T-(ar+afP^)[C,]T-a,r[C,:]T + ar[09j 



dt 



d[C.,] 



kp, [C,_i] + fcdps [C^i] - (fcp. + fcdp.) [C] + [C,]t - [C] 



dt 



-fcf ([B]t - mEjB„.C,]T)"[C,] + XJ 




(30) 




c?[B„iCi]T 
dt 



fcps[BmCi_i]T + fcdps[B,n C+i ] T 



(A;ps + fcdps)[BmCj]T 

g fcf^Kr[B„^C.]T 
K™ -I- [A]™ 



+fcf([B]T 



mEjB„C,]T) 



rn 



(31) 



where the concentration of free KaiA, [A], is given by: 
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N 



[A][CJt 
K. + [A] 




[A]"'[B„,C,]t 
K™ + [A]™ 



[A]t = 



(32) 



The effective (de)phosphorylation rates on the active branch are given by af^ = 
(fcpsKj + fcpf[A])/(K, + [A]) and crfP" = K,fcdps/(Kj + [A]). The dissociation 
constants and are given by = kf^^/kf^^ and K.; — kf^^/kf^. The 
effective flipping rates are given by af^ = fiKi/{Ki + [A]) and af' = hi, where 
fi and bi are the forward and backward flipping rates. We have confirmed that 
for sufficiently large /c^^, k^^ and fc^*, k^^ this set of rate equations gives results 
that are identical to those in Egs. l^TII^ Unless indicated otherwise, the results 
in the main text are obtained by numerically solving Eqs. I29II31I 

3.5 Bifurcation analysis 

We have performed a bifurcation analysis of the full model. To this end, we 
use Eqs. [^Tll^ However, the system of differential equations in Eqs. [^Tll^ 
obeys the conservation law + [ACi] + [Ci] + [B,„Ci] + [AmBmCi]) = [C]t. 

As a consequence, a linear stability analysis would always yield at least one 
eigenvector with eigenvalue zero, which complicates detection of bifurcation 
points. To eliminate the zero eigenvalue associated with the conservation of 
KaiC, we express the concentration of one of the KaiC complexes in terms of 
the concentrations of the other KaiC complexes. We have chosen [AC at] to 
take this role. Thus, [AC at] is not a separate dynamical variable, but is instead 
defined by: 



[AC^] = [C]t-5] [AC2]-E([C,] + [C] + [B„C,] + [A™B™C,]) (33) 



Numerical continuation of the fixed points and limit cycles was performed with 
the software package XPPAUT ^ , which incorporates the numerical continua- 
tion routines from AUTO [S]. 

In Fig. [5] we show the bifurcation diagram of the full model, defined by 
Eqs.HnSni as a function of [A]t, for different values of [B]t; the phase diagram 
is shown in Fig. [51 For very small KaiB concentration, [B]t < 0.05[C]t, the 
system has a single, stable fixed point for all [A]t (Fig. O. For higher KaiB 
concentration, 0.05 < [B]t/[C]t < 0.6, the system is bistable for a range of [A]t 
(see Figs. [5]and[6l): it has two stable steady states and one unstable steady state, 
corresponding to different degrees of KaiC phosphorylation. At the boundaries 
of this bistable region, a stable and unstable fixed point merge via a saddle-node 
bifurcation (Fig. O. We discuss the origin of this bistable regime in more de- 
tail below. For even higher KaiB concentration, 0.6 < [B]t/[C]t < 1-2, one of 
the two stable fixed points, namely that with the lower phosphorylation level, 
becomes unstable for a range of KaiA concentrations. This stable fixed point 
becomes unstable via a supercritical Hopf bifurcation and gives rise to a limit 
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Figure 5: Bifurcation diagram of the full model of Eqs. [2TH251 as a function of 
[A]t, for different values of [B]t- Stable fixed points and unstable fixed points 
are indicated by solid black lines and dashed grey lines, respectively. When the 
fixed points change stability, either saddle-node bifurcations (squares) or Hopf 
bifurcations (circles) occur. Apart from [A]t and [B]t, all other parameters are 
as shown in Table [T] 

cycle. Thus, in this window of KaiA and KaiB concentrations, the system has 
one stable fixed point at high phosphorylation level and a limit cycle. For yet 
larger KaiB concentrations, [B]t > 1.2[C]t, the system has only one unstable 
fixed point surrounded by a limit cycle for a range of [A]t; again, the limit cycle 
appears and disappears at low and high [A]t, respectively, via a supercritical 
Hopf bifurcation. This limit cycle corresponds to the circadian oscillations dis- 
cussed in the main text. Fig. [5] shows that this oscillatory regime with only one 
limit cycle has a lower and an upper bound on the KaiA concentration, but no 
apparent upper limit on the KaiB concentration. In contrast, both the bistable 
regime and the regime in which a limit cycle coexists with a stable fixed point, 
occur only over a fairly narrow range of KaiA and KaiB concentrations. 

In Fig. [3 we examine the properties of the limit cycle that is created at 
the Hopf bifurcation. It is possible to do numerical continuation of the limit 
cycle in the vicinity of the Hopf bifurcation as is shown in Fig. [7l\. This analysis 
shows that the limit cycle is stable and that the bifurcation is thus supercritical. 
Further away from the Hopf bifurcation the numerical continuation algorithm 
fails to converge. However, as shown in Fig. [7J3 and C, by directly solving the 
differential equations I^TH^ we can nonetheless show that the system continues 
to converge to a stable limit cycle. Because of the fact that the algorithm cannot 
continue the limit cycle all the way from one Hopf bifurcation to the other, we 
cannot strictly rule out the possibility that it undergoes further bifurcations. 
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Figure 6: Phase diagram of the fuU model. In the region enclosed by the dashed 
grey lines, the system possesses a stable limit cycle. In the region enclosed by 
the solid black lines, the system has three fixed points, of which two are stable 
in the absence of a limit cycle. Where the two regions overlap, a single stable 
fixed point coexists with a limit cycle. 

Nevertheless, we find by direct integration of the differential equations that both 
the period and amplitude of the limit cycle vary smoothly with [A]x between 
the Hopf bifurcations. 

We now discuss the origin of bistability for intermediate [A]t and [B]t (see 
Fig. ini). Fig. [5] shows four typical time traces of the phosphorylation level 
of KaiC when the system is in the bistable regime. The different time traces 
correspond to different initial conditions. These initial conditions differ in the 
phosphorylation level of KaiC. Indeed, the initial degree of phosphorylation 
largely determines which one of the stable fixed points the system converges to. 
For low initial phosphorylation, the system converge to a steady-state phospho- 
rylation level of — 0.5, while for high initial phosphorylation, it converges 
to a phosphorylation level of = 0.9. These two steady states do not only 
differ in the average phosphorylation level of KaiC, but, importantly, also in the 
concentration of /ree KaiA: for ps — 0.5, [A] is small, while for ps — 0.9, [A] is 
large. 

To understand the origin of the difference between the two steady states, it 
should be realized that: a) KaiB is needed to stabilize inactive KaiC, but its 
concentration is fairly low; b) KaiA is needed to phosphorylate active KaiC, but 
also its concentration is fairly low. 

In the low ps state, most KaiC hexamers have initially a low degree of 
phosphorylation. These hexamers will bind KaiA, which will stimulate their 
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Figure 7: Limit cycle in the full model for [B]t/[C]t = 3. (A) Bifurcation di- 
agram of [C4] in the vicinity of the Hopf bifurcation as obtained by numerical 
continuation of the limit cycle. The stable and unstable fixed points are indi- 
cated by a solid black line and a dashed grey line, respectively. The minimum 
and maximum values of [C4] along the limit cycle are shown as thick black lines. 
The limit cycle is stable, indicating a supercritical Hopf bifurcation. Here, we 
choose to plot [C4] for convenience, and the concentrations of other components 
of the system show similar behavior close to the Hopf bifurcation. (B) and (C) 
Limit cycle for [A]t/[C]t — 1 and [B]t/[C]t = 3, obtained by numerical inte- 
gration of Eqs. [2TII25I (B) Phosphorylation p in time. (C) Phosphorylation p 
versus c, the fraction of KaiC in the active state. 

phosphorylation. However, because [A]t is low, the phosphorylation rate per 
hexamer will be low and counterbalanced by the spontaneous dephosphorylation 
rate. As a consequence, the overall phosphorylation level will be low. 

In the high ps state, most KaiC hexamers go through the approximate cycle 
AC5 — > Cg ^ B2C6 B2C5 AC5. In the high ps state, most KaiC hexamers 
initially have a high degree of phosphorylation. The available KaiA dimers 
will be able to fully phosphorylate these hexamers before they flip towards 
the inactive state. On the inactive branch, these hexamers need KaiB to be 
stabilized. However, because [B]t is low, the inactive hexamers will not be 
stabilized very strongly, and will therefore flip back towards the active state. 
At this point, the concentration of free KaiA is still high, because there are no 
hexamers with a low phosphorylation level, which could bind KaiA. Because 
the concentration of free KaiA is high, the hexamers that have just flipped back 
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Figure 8: Bistability in the full model for [A]t/[C]t = 0.48 and [B]t/[C]t = 0.5. 
The degree of phosphorylation in time is shown for different initial conditions, 
[Ci(0)] — [C]t for i = 1,2,3,4. For low initial phosphorylation, the system 
converges to a steady state at p « 0.5. For high initial phosphorylation the 
system converges to another steady state at p ^ 0.9. 



towards the active state can be rephosphorylated, and the cycle starst again. 
This situation clearly elucidates the important role of KaiB. Without KaiB the 
inactive branch is not stable, and the full allosteric cycle will be short cut. This 
will kill the capacity of the system to generate macroscopic oscillations. 

The above mechanism for generating bistability is similar to the mechanisms 
that have been proposed for generating bistability in the MAPK ^ and the 
CAMKII system [7]. Both in these systems and in the Kai system, a protein can 
be phosphorylated at multiple sites and the limiting amount of cither the kinase, 
as in the Kai or the MAPK system, or the limiting amount of phosphatase, as 
in the CAMKII system, makes the system flip-flop between a state with a high 
phosphorylation level and one with a low phosphorylation level. 

Fig. [S] summarizes the behavior of the full system. The full model has a limit 
cycle for a broad range of concentrations. Although the range of [A]t for which 
oscillations are observed decreases slightly with increasing [B]t, we found no 
indication that oscillations cease for higher [B]t. This is in agreement with the 
passive role played by KaiB in stabilizing the inactive branch and sequestering 
KaiA. The region for which bistability occurs is much smaller. Furthermore, 
as we have discussed above, the occurrence of the bistable regime does depend 
upon details of the model, such as the extent to which KaiB sequesters KaiA 
and stabilizes the inactive state when bound to KaiC. Further experiments will 
be needed to determine whether bistability really occurs in the Kai system. 
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